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We present an analytic method to determine spectral properties of the covariance matrices con- 
structed of correlated Wishart random matrices. The method gives, in the limit of large matrices, 
exact analytic relations between the spectral moments and the eigenvalue densities of the covariance 
matrices and their estimators. The results can be used in practice to extract information about the 
genuine correlations from the given experimental realization of random matrices. 

Wishart random matrices play an important role in the multivariate statistical analysis 0, Q • They are useful in 
some problems of fundamental physics Q , communication and information theory 0, Q , internet trading Q and 
quantitative finance |a l9l flfj| . 

A Wishart ensemble of correlated random matrices is defined by a Gaussian probability measure: 



P{X)DX = A^exp 
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where X = (Xi a ) is a real rectangular matrix of dimension N xT. It has two types of indices: an A-type index runnig 
over the set i — 1, . . . , A and a T-type index over a — 1, . . . , T. Throughout the paper the A- type indices will be 
denoted by Latin letters and the T-type by Greek ones. X r denotes the transpose of X. The matrices C = (CV,) and 
A = (A a p) are symmetric square matrices of dimensions A x A and T x T, respectively. They are positive definite. 
A" is a normalization constant: 

N= (2^)^f(detC)*(detA)Tr ; (2) 

chosen to have J P(X)DX = 1. 

Let Q(X) be a quantity depending on X. The average of Q over the random matrix ensemble is defined as: 

(Q) = j Q(X)P(X)DX. (3) 

In particular, the two-point correlation function is: 

(X ia Xjp) = CijA a p, (4) 

as directly follows from the Gaussian integration. In this paper we are interested in the spectral behaviour, the 
eigenvalue distribution and the spectral moments of the following random matrices: 

c=ixX^ , a=ix T X. (5) 

These matrices can be used as estimators of the correlation matrices C and A if some realizations of random matrices 
X are given. We will refer to c and a as to covariance matrices or statistically dressed correlation matrices. We will 
present an analytic method to determine the eigenvalue distribution and the spectral moments of c and a in the limit 
of large matrix size. Another method of calculating the eigenvalue density of correlated Wishart matrices has been 
recently discussed in [Tlj . 

In parallel to Q one can define a Wishart ensemble of correlated complex matrices: 

V(X)DX = Af' 1 exp [-TrXtC^XA- 1 ] JJdA^dA™. (6) 
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The matrices C and A are now Hermitean and positive definite. denotes the Hermitean conjugate of X. The 
normalization constant is now AT = (7r) ArT (det C) T (det A) N . In the analysis of the complex ensemble the estimators 
of the correlation matrices are replaced correspondingly by 

c = ixxt , a=lxtX. (7) 

Notice that the factor one half in front of the trace in the measure for real matrices Q is dropped in JJJJ. With this 
choice of the measure the two-point correlations take a similar form as for real matrices Q : 

{Xi a X*p) — CijA a p . (8) 

The star stands for the complex conjugation. Additionally we also have: 

(x ia x jl3 ) = (x* a x; p ) = o . (9) 

As a consequence, as we shall discuss towards the end of the paper, the matrices c and a © in the real ensemble have 
an identical large N behaviour as the corresponding matrices Q) in the complex ensemble. Since we are interested 
here only in the large N behaviour it is sufficient to consider one of the two ensembles and draw conclusions for the 
other. We will focus the presentation on the ensemble of real matrices. 

An example of a problem which can be formulated in terms of Wishart random matrices Q is the following. 
Imagine that we probe a statistical system of N correlated degrees of freedom by doing T measurements. We store 
the measured values of the i-th degree of freedom in the a-th measurement in a rectangular matrix X = (Xi a ). The 
degrees of freedom as well as the measurements may be correlated. This is expressed by the equation (@}, which tells 
us that the covariance matrix for the correlations between degrees of freedom in the system is C — {Cij ) and for the 
(auto)correlation between measurements is A = (A a p). Note that in general the correlations between Xi a and Xjp 
may have a more complicated form: (Xi a Xjp) = Cia^p, where the matrix C has double indices. Such a situation 
takes place if the autocorrelations are different for various degrees of freedom. We shall not discuss this case here. 
Moreover, we shall assume that only Gaussian effects are important for the studied system. 

A perfect example of the situation described above is the problem of optimal portfolio assessment - one of the 
fundamental problems of quantitative finance. The portfolio assessment is based on the knowledge of the covariance 
matrix C for stocks returns [l^ . In practice, the covariance matrix is estimated from the historical data which are 
stored in a rectangular matrix representing T historical values of N stocks. Fluctuations of returns are well described 
by the Gaussian ensemble . The estimator of the covariance matrix is given by (0 ■ Another problem of modern 
financial analysis which can be directly cast into the form is the problem of taste matching 0. This problem is 
encountered for instance in the large-scale internet trading. 

It is worth mentioning that the random matrix framework may also be used in a statistical description of data 
generated in Monte Carlo simulations for a system with many degrees of freedom, in particular of data concerning the 
correlation functions. One frequently encounters such a problem in Monte-Carlo simulations of lattice field theory, 
where the field is represented by correlated numbers distributed on a lattice. Usually, one is forced to use a dynamical 
Monte Carlo algorithm to sample such a system. The basic idea standing behind a dynamical algorithm is to generate 
a Markov chain - a sort of a random walk - in the space of configurations. The degrees of freedom on the lattice as 
well as the successive configurations are usually correlated. Outside a critical region no long range correlations are 
observed and the fluctuations can be treated as Gaussian. 

Complex random matrices are useful for instance in telecommunication or information theory 

Let us come back to the ensemble of real matrices . As we mentioned the matrices (0 can be treated as estimators 
of the correlation matrices C and A. Indeed, from equation (J3J we see that: 



(en) = ( \ E x *° x ^ ) = M ^ c v > ( 10 ) 

{a a p) = (^y^ y X ia Xip\ = MciAap, (11) 

where Mci = ^-Tr C and Mai = ^Tr A. This notation will be explained later. The last equation tells us that 
measuring the average of c over the ensemble Q we obtain the matrix C up to a constant. In other words having 
a realization of random matrices X we can use (jSJ to estimate C. Similarly, we can use a to estimate A. Notice 
that the measure Q is invariant under the transformation C — > bC and A — * 6 _1 A, where b is an arbitrary positive 
real number. In particular (cy) and {a a p) are independent of the rescaling factor b. This independence is ensured by 



the presence of the factors Mai and Mci in equations i|10lll|) . In practical calculations, if Tr A and TrC are not 
specified, one can remove the redundancy with respect to the rescaling by b, setting yTr A = -^TrC. In this case 
the constants Mai = Mci can be determined from the data by evaluating the traces of c or a: 



Moi - Mai = \J j^Tr (c) = ^ ^Tr (a) . 

While considering the covariance matrices for the Wishart ensemble we can formulate two reciprocal problems, 
which we shall call direct and inverse problem. In the direct problem we want to learn as much as possible about 
the probability distribution of the estimators c and a JSJ assuming that the matrices C and A are given. In particular, 
we want to calculate the eigenvalue density functions: 

Pa(A) = ^Xj£(A-A a )\ , 

where and A Q are eigenvalues of c and a, respectively. The determination of the eigenvalue density functions is 
equivalent to the determination of all their spectral moments: 



m ck = J d\p c {\) \ k = ^Trc fc 
m ak = J d\ Pa {\) X k = ^Tra fc 



The moments m c k,m a k are related to each other: 

TOafc = r l - k m ck , (12) 

where r = N/T, as follows from the cyclicity of the trace: 

Tr XX 1 ■ • ■ XX 1 = Tr X r X ■ • • X T X . 

In the inverse problem we want to learn as much as possible about the genuine correlations in the system, which are 
given by C and A, using a measured sample of random matrices X. We can do this by computing the estimators c 
and a J^J) and relating them to matrices C, A. In particular we would like to estimate the eigenvalue distributions 
and the moments of C, A: 



I 1 - 



Mr.k = 4zTrC k = -VA l 



i=\ 
T 

M Ak = iTr A fc = 1 A* , 

a=l 

where Aj and A Q are eigenvalues of C and A, respectively. The inverse problem is very important for practical 
applications, since in practice it is very common to reconstruct the properties of the underlying system from the 
experimental data. 

In the analysis of the spectral properties of the matrices a and c it is convenient to apply the Green's function 
technique. One can define Green's functions for the correlation matrix C and its statistically fluctuating counterpart 
c: 

Gc(*)= ; 1 r ■ (13) 

^) = (-^—), (14) 
\zl N - c 



and correspondingly Ga(z) and g a (z) for A and a. The symbol ljy stands for the N x N identity matrix. A 
corresponding symbol It appears in the definition of Ga(z) and g a (z)- The Green's functions are related to the 



generating functions for the moments: 



z fe AT 
fc=i 



.(*) = £^ = i^( *&(*))-!. (16) 



z fc AT" 
fc=i 



or inversely: 



_L TrGc(2) = i±j^<£> 



-Trg c z = Eii . 17 

TV z 

The analogous relations exist for Ga(z) and g a (z)- The Green's functions can be used for finding the densities of 
eigenvalues: 

Pc(A) = -ilml T r gc (A + i + ) = _I Im i±T^±!^l , (18) 

7T TV 7T A + l(J + 

and similarly for p a (A). The eigenvalue densities p a (A) and p c (A) are not independent. As follows from (|12|l the 
corresponding generating functions (|16J) fulfill the equation: 

w a (z) = rm c (rz) . (19) 

Combining the last equation with (|17|) we obtain: 

^Tr g a (z) = r 2 ^Tr g c (rz) + ^ . (20) 

Applying now Ijl8(l we have: 

Pa(A) = r 2 p c (r\) + (l-r)5(\) . (21) 

The meaning of the last term on the right hand side of this equation is that there are T — TV zero modes in the matrix 
a if T > TV. The zero modes disappear when TV = T. Moving the term containing the delta function to the other side 
of equation, dividing both sides of the equation by r 2 and substituting the parameter r by s = T/N = 1/rwe obtain: 

Pc(A) = s 2 p a (sX) + (l- s )S(X) . (22) 

Therefore, for T < TV the zero modes appear in the spectrum /? C (A). In this case it is more convenient to use the 
parameter s = 1/r instead of r. The zero modes appear in the eigenvalue distribution of either a or c. The two 
equations l|21|l and 12211 are dual to each other. For r = s = 1 they are identical. Because of the duality it is sufficient 
to solve the problem for r < 1. We will present a solution for the limit r = N/T = const < 1 and N —> oo neglecting 
effects of the order 1/TV. 

Using a diagrammatic method 0, Q, 0] one can write down a closed set of equations for the Green's function 
gc(«) G3: 

8c{z) = zl N - E c (z) ' g * c(z)= Tl T -£» c (z) ' 

(23) 

£ c (z) = CTr(Ag„ c (z)) , S, c (z) = ATr (C g c (z)) . 

The set contains four equation for four unknown matrices including g c (z) which we want to calculate, and three 
auxiliary ones: g* c (-z), S c (z), E* c (z) (see Appendix 1). Each of them can be interpreted in terms of a generating 
function for appropriately weighted diag rams with two external lines: g c (-z), g* c (z) for all diagrams and S c (z), S* c (z) 
for one- line-irreducible diagrams jl3Ul4ill!T| (see Appendix 1). In the limit TV — > oo the weights of non-planar diagrams 
vanish at least as 1/TV. Thus in this limit only planar diagrams give a contribution to the Green's function. Therefore 
the large TV limit is alternatively called the planar limit. The diagrammatic equations l|23|) hold only in this limit. An 



analogous set of equations can be written for the Green's function g a (^). The equations are identical to those of 1231) 
if one exchanges a < — > c, A < — > C and T < — > N. The two sets can be solved independently of each other. However, 
as follows from the duality it is sufficient to solve only one of them and deduce the solution of the other. 

The equations l|23[l can be solved for g c (z) by a successive elimination of g* c (z), S c (z) and S* c (z). However, the 
resulting equation is very entangled |T^: 



gcM ^ ; l,-CT t lTr _ ATt(CgcM) j , (24, 

and cannot be easily used in practical calculations of the moments m c fc or spectral density p c (A). 

Another way of solving the equations ()23|1 was proposed in [15J] . It relies on introducing a new complex variable Z 
conjugate to z which is defined by the equation: 

m c (z) = M C (Z) . (25) 

At the first glance this equation looks useless because it refers to an unknown function m c (z) which we actually want 
to determine. Quite contrary to this, as we shall see, the introduction of the conjugate variable Z allows us to write 
down a closed functional equation for m c (z). First, let us illustrate how the method works for A = It In this 
case, the elimination of the auxiliary functions 12: il) leads to 

(26) 



1 + rm c (z) 
or equivalently to 

z = Z ■ (1 + rAl c (Z)) . (27) 

Suppose we solve the direct problem. In this case we know the matrix C and hence also the generating function 
Mc{Z). Inserting l|26|) to (|25|l we obtain a closed compact functional relation for m c (z): 



If Mc(z) has a simple form, one can solve the equation for m c (z) analytically |l5j |. In general, one can write a 
numerical program to calculate the eigenvalue density p c (X) from the last equation. In case of solving the inverse 
problem, we assume that we can determine moments m c k from the data and hence that we can approximate the 
generating function m c (z). Then we can insert l(27|l to (|25J) and obtain a functional equation for Mc(Z): 

M c {Z)=m c (Z-(l + rM c (Z))). (29) 

The problem is solved in principle. However, in practical terms the inverse problem is much more difficult, because one 
cannot compute all experimental moments m c k with an arbitrary accuracy, unless one has an infinitely long series of 
measurements. But one never has. In practice one can estimate only a few lower moments m c k with a good accuracy. 
Because of this practical limitation one cannot entirely solve the inverse problem. However, as we discussed in [T3| 
the inverse problem can be partially solved even in specific practical applications using a moments method. Let us 
sketch this method below. 

We can gain some insight into the spectral properties of the correlation matrix C by determining the relation 
between the moments Met and m c k- Expanding the functions Mc(z) and m c (z) in (|28|l in 1/z using (|15ll(j|) and 
comparing the coefficients at 1/ z k we obtain: 

m cl = Mci 

m c2 — M C 2 + tMq X 

m c3 = M C3 + 3rMciM C2 + r 2 M c 1 (30) 

m c4 = M C4 + 2r (M C2 + 2M C1 M C3 ) + 6r 2 M cl Af C2 + r 3 M^ 



We can also invert the equations for Met. The result of inversion gives a set of equations which can be directly 
obtained from the 1/Z expansion of the functions in the equation (|29|l which is the inverse transform of (|28|l . We can 
also determine the corresponding relations for 'negative' moments m c k and Mcfc, that is for k < 0, or determine the 
spectral density p c (X) [15J. Using a computer tool for symbolic calculations one can easily write a program which 
successively generates the relations between spectral moments (I3U|) from the equation 1)28(1. 



The calculations get more complicated in the general case when both C and A are arbitrary. The guiding principle 
is the same, though. We introduce the conjugate variable Z ((25)) and, using it, write down the solution of the equations 
(j2SJ- In the direct problem we assume that the generating functions Mc{Z) and M A (Z) are known. We will show 
that in this case the solution of (|23[l takes a form of an explicit equation for z = z(Z), where the function z(Z) 
depends on the functions Mc(Z) and M A (Z). Inserting this solution back to (|23[) we eventually obtain a functional 
equation m c (z(Z)) = Mc(Z) from which we can extract the function m c (z). 

The solution of (|23|) takes the form: 

Z r^Jl-A a rfM c (Z)' [ ' 

where A Q are eigenvalues of A. It can be rewritten as: 



rMc ^ = M H^cMJ • (32) 

and can be formally solved for z: 

z = Z rM c (Z) M A 1 (rM c (Z) ) , (33) 

where M A is the inverse function of M A . Thus we have obtained an explicit equation for z = z(Z) in terms of the 
known functions M A and Mc- One can easily check that for A = It the last equation reduces to (|27|l . In this case 
M A (z) = l/(z - 1), M A \z) = 1 + l/z. 

Combining the equation for z = z(Z) given by l|33|) with (|25|l we arrive at a closed equation for the generating 
function m c (Z). It can be used for example to calculate the moments to^'s (see Appendix 2). The calculations yield 
a set of equations expressing m c fc's in terms of the bare moments M A k and Mck- 

m cl = M C1 M A1 

m c2 = M C2 M 2 Al + rM^M A2 

m c3 = M C3 M A1 + 3rM C iM C2 M A1 M A2 + r 2 M^M A3 (34) 
m c4 - M C4 M A1 + 2r (M£ 2 + 2M C iM C3 ) M 2 Al M A2 + 2r 2 M 2 cl M C2 (M A2 + 2M A iM A3 ) + r 3 M^M Ai 

The equations reduce to the form (|30() for A = 1^. 

Using the relations 1|12J) we can also determine the moments of the matrix a. It is more convenient to write them 
using the variable s — r~ x - the dual counterpart of r - instead of r itself: 

TOal = M A1 M C1 

m a2 = M A2 M 2 C1 + sM AX M C 2 

m a3 = Af A3 A^ci + 3sM A1 M A2 M cl M C 2 + s 2 M A1 M C3 (35) 
m a4 = M A4 M^ + 2s (M| 2 + 2M A iM A3 ) M^Mcz + 2s 2 M A1 M A2 (M^ 2 + 2M C1 M C3 ) + s 5 M A1 M C i 

The equations are completely symmetric to Ij34(l with respect to the change r < — ► s (which amounts to N < — > T), 
and c < — > a, C < — > A. Using this method one can obtain equations Ij34(l and 1)35(1 to an arbitrary order. 

The above relations are useful for computing the dressed moments w a fc,?7icfc for given matrices A, C or inversely, 
the genuine moments M A k, Met from the experimental data. As mentioned, the spectral moments give us in principle 
full information about the eigenvalue distribution. In practice the reconstruction of the eigenvalue density may be 
difficult, because to do it we would need to know all moments with a very good precision. Usually, in practical 
applications one can accurately evaluate only a few lower moments. 

In some special cases if we can make some extra assumptions about the form of the matrices C or A we can 
improve significantly the reconstruction of the eigenvalue density. In the previous work [l5l | we have analysed the case 
of A = It and of the matrix C which had only a few distinct eigenvalues. In this case the Green's function g c (^) is 
given by an algebraic equation of the order which is equal to the number of distinct eigenvalues. It can be analytically 
solved when this number is less or equal four. If it is larger the problem can be handled numerically. The duality tells 
us that the solution also holds when we change the roles of A and C. 

Below we will discuss the case of exponential autocorrelations. Exponential correlations are encountered in many 
situations. The general solution, which we have discussed so far, simplifies in this case to a more compact relation 



for the Green's function, which allows us to find analytically an approximate form of the eigenvalue density of the 
random matrices c and a. The approximation becomes exact in the large N limit. We consider purely exponential 
autocorrelations given by the autocorrelation matrix: 



A a p = exp [-\a-f3\/t] , 

where t controls the range of autocorrelations. The inverse of the matrix A reads: 

ex, —1, 
-1, 2ch, -1, 



(36) 



2sh 



-1, 2ch, -1 
— 1, ex 



(37) 



We have introduced here a shorthand notation ex — exp(l/t), ch — cosh(l/i) and sh = sinh(l/£). The spectrum of 
this matrix can be approximated by the spectrum of a matrix M: 



M 



2sh 



2ch, -1, 
-1, 2ch, -1. 



-1, 2ch, -1 
-1, 2ch 



(38) 



whose eigenvalues can be found analytically: 



fi a = (ch + cos(ira/T))/ sh. 

The corresponding eigenvectors are given by the Fourier modes. The matrix 2sh ■ M can be viewed as of a sum: 
(2ch — 2) + At, of a unity matrix multiplied by a constant and a discretized one-dimensional Laplacian At for a 
cyclic chain of length T. The matrix A -1 can be obtained from M by adding to it a perturbation P: A -1 = M + P, 
where P has only four non-vanishing elements: Pn = Ptt = ex -1 and Pit = Pti = 1- The first order corrections 
to the eigenvalues of A -1 , which stem from the perturbation P, behave as 1/T. The perturbation P can be viewed 
as a change of a boundary condition of the Laplacian. As usual, boundary conditions affect mostly the longest 
(small momentum) modes. Indeed, a careful analysis shows that the two diagonal terms of the perturbation matrix, 
Pn = Ptt, introduce a constant correction independent of T of the lowest eigenvalues which does not vanish when 
T goes to infinity. However, since the differences between unperturbed eigenvalues of M and the corresponding 
perturbed eigenvalues of A -1 disappear for all other eigenvalues, we expect that for t <^_T the spectral properties of 
A can be well approximated by the eigenvalues of M _1 : 



fi a ch + cos(na/T) 
In this limit we can also approximate the sum Q31[l by an integral: 



(39) 



z sh 



Z ir J (ch — sh ■ F) + cost 



(40) 



Where t = -Kct/T and the symbol F = ^rm c (z) is introduced for brevity. Note that in the definition of F we have 
replaced Mc(Z), which would rather be dictated by (fTT3~f> . by m c (z). This change is legitimate due to (|2*5|l . The 
integral IpEtH) can be done: 



sh 



Z ^(ch-shF) 2 - 1 



(41) 



Setting back F = —rm c (z) we eventually obtain: 



Z = z 



ch ■ rm c (z) — yl sh 2 + r 2 m\(z) 
sh ■ (r 2 m 2 (z) — 1) 



(42) 



t = 1 




t = 5 






t 


= 10 




T 


M A 2 


Mas 


M A 4 


T 


M A2 


Mas 


M A4 


T 


M A2 


Mas 


Ma4 


20 


1.29493 


2.01479 


3.48620 




20 


4.44996 


28.6455 


204.107 




20 


7.58726 


79.6134 


891.336 


50 


1.30579 


2.05757 


3.60838 




50 


4.81980 


34.2544 


271.932 




50 


9.03668 


120.509 


1784.56 


100 


1.30941 


2.07183 


3.64911 




100 


4.94314 


36.1292 


294.733 




100 


9.53497 


135.501 


2146.93 


200 


1.31123 


2.07896 


3.66947 




200 


5.00482 


37.0666 


306.133 




200 


9.78414 


143.001 


2328.48 


500 


1.31231 


2.08324 


3.68169 




500 


5.04182 


37.6290 


312.973 




500 


9.93364 


147.501 


2437.40 


oo 


1.31304 


2.08609 


3.68983 




oo 


5.06649 


38.0040 


317.534 




oo 


10.0333 


150.501 


2510.02 



TABLE I: The moments Ma2, Ma.3 and Ma4 of the matrix A l|36fl for three values of the autocorrelation length t = 1,5, 10 
calculated numerically for finite size T = 20, . . . , 500 and by the analytic formula 1451 which corresponds to T — oo. The finite 
size values approach the values given by 14511 as T/t tends to infinity. 



This is an explicit equation for Z = Z(z) which can be now inserted into m c (z) — Mc(Z) giving us a compact relation 
for m c (z) in the presence of the exponential autocorrelations 1|36JI : 



/ ch ■ rm c {z) - y/ sh 2 + r 2 m 2 c {z) \ 
mc(z) = MC [* sh.(r*ml{z)-l) ) ■ (43) 

In the limit t — > 0, the parameters sh = sinh(l/i), ch = cosh(l/i) and ex — exp(l/t) increase to infinity and 
ch/ex w sh/ex w 1. As a consequence, the form of equation 14211 simplifies to (|26|l . which corresponds to the 
case without autocorrelations, as expected. Using the equation (|43|l we can recursively generate equations for the 
consecuitive moments: 



m c i 


= Mci 




= Mc2 


m c3 


= M C3 


m c4 


= M C4 



„ 2 fl/f3 



3,9 1 

- cth 2 - - 
2 2 



(44) 



where cth = ch/sh = coth(l/£). The coefficients on the right hand side, which depend on cth, can be directly 
expressed in terms of the moments M^k of the matrix A. Approximating again a sum by an integral in the large T 
limit we can write: 



M 



Afe 



rp 7T 

T ^ a IT J 



(It 



The integrals can be calculated yielding: 



T ^ a it J (ch + cosr) 

a. — 1 n 



M A i =1 
Ma 2 = cth 



M A3 = \cth 2 - \ (45) 
M A 4 = \cth z - \cth 



We see that if we insert these coefficients into the equations (I34JI we obtain (|44(l . This is a consistency check for the 
approximation which we use here. The quality of this approximation can also be checked by comparing the moments 
M\k of the matrix (|36(l for finite T with the result l|45(l which corresponds to T = oo. We expect that for t <C T the 
numerical values shall approach the result I45|l . The results of this comparison confirm our expectations (see table 1). 
Thus we see that the formula 144(1 for m c (z) in the presence of the exponential autocorrelations becomes exact in the 
limit T — > oo. This formula allows us to compute the eigenvalue distribution of the random matrices c and a (J3J). Let 
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FIG. 1: The density of eigenvalues p c (A) for exponential matrix A, C 
times t: t = (solid line), t — 2 (dashed line) and t = 5 (dotted line). 



ljv,r = 0.2 and for three different autocorrelation 



us illustrate this on the simplest example of the system which has no correlations: C = ljv and Mc(Z) = l/(Z — 1) 
In this case the equation i|43|l for m c (z) takes the form: 



cth ■ x 2 (x 2 - l)(x + r) - x^Jl + x 2 /sh 2 = 



where we have used the notation x(z) = rm c (z). For t — * (A 



-(x + l)(x + r) 



It) this equation reduces to: 
. 



which has a solution: 



x(z) = - [-1 -r- iy/{n+ - z){z - /!_) + z j 
where (j,± = (1 ± \fr) 2 ', which leads to the well-known result for the uncorrelated Wishart ensemble: 

Pc(A) 



1 v/( M+ -A)(A-/i_) 



2irr A 

For t > it is still possible to find analytically a solution of Q46|). Let us rewrite (|46|) as a polynomial equation: 



(46) 



(47) 



(48) 



(49) 



z 2 x 2 {l + x 2 /sh 2 ) + (cthzx 2 - (r + x)(x 2 - 1))" 







It has two trivial solutions x = ±1. Dividing out the polynomial (x — l)(x + 1) we get: 

x 4 + 2x 3 (r - cthz) + x 2 (-l + r 2 - 2cthrz + z 2 ) - 2rx - r 2 = 







(50) 



This is a quartic equation which can be solved analytically by the Ferrari method. We will not present the formal 
solution which is neither transparent nor informative. Instead, we show in Fig. ^ the eigenvalue density functions 
p c (A), for different t, resulting from this solution. The lower part of the distribution approaches zero when t increases, 
but zero modes do not appear in the distribution as long as r < 1. The formula l|43(l applies to any correlation matrix 
C but in the general case one has to use a numerical procedure to calculate from it the density function. 

Let us stop here the presentation of results for the ensemble of real matrices. As we mentioned all results in the 
large N limit hold also in the ensemble of complex matrices if the covariance matrices (JSJ) are replaced by Q. The 
reason why it is so is related to the fact that the moments of c = yXX^ in the ensemble of complex matrices © are 
equal to the moments of c = ^XX T in the real ensemble up to a 1/N corrections which disappear in the large N 
limit: 




complex 




0(1/N) 



(51) 



real 



Let us illustrate this by explicit calculations of the second moment. Using the Wick's theorem for Gaussian integrals 
and the equation J3J for the two-point correlation function, we have: 



N 



1 




;XX T J J - (XiaXjaXjpXip) - 

{(Xi a Xj a ) (XjfiXip) + {XiaXifj) (Xj a Xjp) + {X ia Xjfj) (Xj a Xip)} 




M C iM\ x + tMI x M A 2 + —M C2 M A2 . 



The corresponding calculations for the complex ensemble read: 




N 



1 



( 



-XX^ ^ - (X la X* a X 3l3 X* ) - 

2 {(XiaX* a ) (XjpXip) + (X ia X* (j ) (XjaXjp) + (XiaXjp) (X* a X*p)} 



= M C2 M A1 + rMl 1 M A2 . 



The difference between the two calculations appears in the third term which in the real ensemble gives a contribution 
of the order 1/iV while in the complex ensemble disappears by virtue of 10. We recognize that Mc 2 M A1 + rM^, 1 M A2 
which are the leading terms in the l/N expansion are identical as in the second equation in the set (|34[l . Generally 
one can show that the leading contributions which correspond to the planar diagrams in the expansion of gc(z) are 
identical for both ensembles. Non-planar diagrams are different but they contribute in the subleading orders: it turns 
out that in the diagrammatic expansion of the Green's function g c {z) for the complex matrix ensemble ©, which 
would be a counterpart of 1521) in the appendix, all diagrams which contain a double arc with dashed and solid line 
crossed are identically equal zero since such an arc corresponds to the propagator {XiaXjp) or (X* a X*p) . A crossing 
of two arcs is however allowed and leads to a factor l/N 2 . 

To summarize: in the paper we have considered an Wishart ensemble of correlated random matrices. We have 
obtained in the limit of large matrices a closed set of equations relating the Green's function or equivalently the 
moments' generating functions m c (z) and m a (z) for statistically dressed correlations to the generating functions for 
genuine correlation matrices Mc(z) and M A (z). The equations in the large N limit are the same for the ensemble of 
real and complex matrices. Using these equations we can write down exact relations between genuine and experimental 
spectral moments of correlation functions of an arbitrary order. The relations can be used in practical problems to 
learn about correlations in the studied system from the experimental samples. In the case of exponential correlations 
we have also found an explicit form the spectral density function of the covariance matrix. A natural generalization 
of the work presented here is to consider a more general type of time correlations than purely exponential (|36|) . If 
the correlations are of the form which depends on the time difference A a ^p = A(\a — (3\) and if they are short-ranged 
then one can apply Fourier transform to determine in the large T limit an approximate spectrum of the matrix A 
and approximate values of its spectral moments. Another interesting issue which can be addressed in the future is 
the determination of the probability distribution for individual elements of the covariance matrices c and a similarly 
as it was done for the uncorrelated case [T3 |. 
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Appendix 1 



For completeness we recall here the graphical representation of the Green's function. The details of the diagrammatic 
method can be found in 0,0,0]. The Green's function g c (z) can be represented as a sum over diagramms: 




+ 




+ ■■• (52) 

where the iV-type and T-type indices of X are denoted by filled and empty circles, respectively. The matrix X is 
denoted by an ordered pair of neighbouring filled and empty circles, while X T is drawn as an pair of such circles in the 
reverse order. A horizontal solid line stands for ljf/z, a dashed line for 1t/T, a solid arc for C and a dashed arc for 
A. The two point function @ is drawn as a double arc. Matrices on a line arc multiplied in the order of appearance 
on this line. If a line is closed, the trace is taken. 

In the thermodynamical limit only planar diagrams give contribution to g c . In particular the last term in Ij52(l 
vanishes. The Green's function g c *(z) is represented by an identical set of diagrams with dashed and solid lines 
exchanged. It is convenient to introduce one-line irreducible diagrams and corresponding generating functions S c and 
X* c . The Green's functions can be expressed in terms of S c and X* c as follows: 



(53) 



o- 



- -o 



■ o o--o o--o o--o o--o o--o 



(54) 



In the planar limit there are two additional equations which relate the sums over one-line irreducible diagrams to the 
Green's functions: 

(55) 

Analogous diagrammatic equations can be written for g a with the only difference that the solid line shall denote the 
propagator 1t/z and the dashed line l^/N. 





Appendix 2 



We use the equations (|33|l and (|25(l to determine the relations . As for the case A = It we shall do this using 
1/z expansion. The function M A (Z) is given by the series: 



M A {Z) 



M A i , M A2 M. 



A3 



z z 2 z* 

r-1 



Let us determine the expansion for the inverse function M A as a series around zero: 

M A \y) = Mauj- 1 (1 + f ny + f W + ...). 
The coefficients of the series can be directly calculated from the condition: 

y = M A (M A \y)) , 

which gives us: 



M. 



Mi 



A2 



(Mai) 



M2 



M A3 M A1 - (M A2 ) 2 
(Mai) 4 



(56) 

(57) 
(58) 

(59) 



The equation (|33|l takes the form: 



z = M A1 Z (1 + tarMc(Z) + n 2 r 2 M%(Z) + . . . ) , 



(60) 



or if written for l/z: 



1 



1 1 



(1 - furMciZ) + (nl - ^)r 2 M 2 c {Z) + ...), 



(61) 



z 



Mai Z 



where: 



M C (Z) 



Max Mc2 Mcs 

z z 2 z 3 



+ 



(62) 



Thus we have expressed l/z as a series of 1/Z. Inserting this series to the equation (125(1 and comparing coefficients 
at 1/Z k we eventually obtain ((34(1 . 
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